ANOVA conventions may not work as you expected

Preface

This is the analysis codes used in the project of “ANOVA conventions may not work as you expected.”

1 Reviewing ANOVA conventions

This section is the analysis of reviewing ANOVA prevalence and complexity, as well as the prevalence of ANOVA conventions.

Package loading and general settings:

# load packages
library(tidyverse)
library(vcd)
library(ggpubr)
library(RColorBrewer)

theme_set(papaja::theme_apa())

# set global chunk options, put figures into folder
knitr::opts_chunk$set(
  echo = TRUE, warning = FALSE, message = FALSE, 
  include = TRUE, tidy = FALSE, archive = FALSE,
  size = "big", fig.width=8, fig.asp=0.7, 
  width = 1800,
  emmeans=list(msg.interaction=FALSE))

source(here::here("R", "funcs.R"))
source(here::here("R", "func_kappa.R"))

# colors used in figures
two_colors <- c("#D55E00", "#56B4E9")
custom_clr <- c("#7F7F7F", 
                "#E64B35", "#56B4E9", "#DF8F44", 
                "#008B45", "#925E9F", "#415384", 
                "#000000")

set.seed(2022)

Only load paper IDs and the journals for later use:

# article ID and journal
df_publ <- readxl::read_excel(here::here("data", "CodingSheet.xlsx"),
                              sheet="PaperMeta",
                              range="A1:G750") %>%
  transmute(`Paper ID`, journal = `Publication Title`) %>% 
  mutate(journal = factor(journal),
         journal = fct_recode(journal,
                              JEPG="Journal of Experimental Psychology: General",
                              `Psych Sci`="Psychological Science",
                              `J Abn Psych`="Journal of Abnormal Psychology",
                              JCCP="Journal of Consulting and Clinical Psychology",
                              JESP="Journal of Experimental Social Psychology",
                              JPSP="Journal of Personality and Social Psychology"),
         journal = fct_relevel(journal, "JEPG", "Psych Sci", "J Abn Psych", "JCCP", "JESP", "JPSP")) 

1.1 Intra-rater consistency

Estimate the intra-rater consistency for the information on whether the study employed ANOVA, experimental designs and the ANOVA prevalence conventions inspected by Y.J. and H.J. separately.

1.1.1 Research article and ANOVA

Load the data of title pre-screening and research articles:

# "Title_Prescreening" Sheet
df_prescreen <- readxl::read_excel(here::here("data", "CodingSheet.xlsx"),
                                   sheet="Title_Prescreening",
                                   range="A1:C750") %>%
  left_join(df_publ, by="Paper ID") %>%
  select(-Title)

# "ResearchPaper" Sheet
df_mainscreen <- readxl::read_excel(here::here("data", "CodingSheet.xlsx"),
                                    sheet="ResearchPaper",
                                    range="A1:L750") %>%
  left_join(df_publ, by="Paper ID") %>%
  select(-Title)

Calculate Cohen’s Kappa of “whether it is a research article”:

# whether it is a research article (unweighted)
df_res_K <- df_mainscreen %>%
  select(`Research article_YJ`, `Research article_HJ`) %>%
  group_by(`Research article_YJ`, `Research article_HJ`) %>%
  summarize(count = n(), .groups = "drop") %>%
  pivot_wider(names_from=`Research article_HJ`, values_from = count, values_fill = 0)

Ka_res <- df_res_K %>%
  select(-`Research article_YJ`) %>%
  as.matrix() %>%
  Kappa()

print(Ka_res, CI=TRUE) # Kappa for research articles

##             value     ASE     z Pr(>|z|)  lower upper
## Unweighted 0.9713 0.01654 58.72        0 0.9389 1.004
## Weighted   0.9713 0.01654 58.72        0 0.9389 1.004

Calculate Cohen’s Kappa of “whether ANOVA was used in the research articles”:

# whether ANOVA was used (unweighted)
df_aov_K <- df_mainscreen %>% 
  filter(`Research article`==1) %>% # only among research articles
  select(`Use ANOVA_YJ`, `Use ANOVA_HJ`) %>% 
  group_by(`Use ANOVA_YJ`, `Use ANOVA_HJ`) %>% 
  summarize(count=n(), .groups="drop") %>% 
  pivot_wider(names_from=`Use ANOVA_HJ`, values_from = count, values_fill = 0) 

Ka_aov <- df_aov_K %>% 
  select(-`Use ANOVA_YJ`) %>% 
  as.matrix() %>% 
  Kappa()

print(Ka_aov, CI=TRUE) # Kappa for using ANOVA

##             value     ASE    z Pr(>|z|) lower upper
## Unweighted 0.8685 0.01965 44.2        0  0.83 0.907
## Weighted   0.8685 0.01965 44.2        0  0.83 0.907

1.1.2 The first 120 articles

When trying to estimate the consistency of analysis information for the first 20 articles in each journal, we realize that we may not identify the same number of ANOVAs for each experiment/article. Therefore, we add some un-preregistered analysis steps. Specifically, before evaluating the consistency for the ANOVA information, we count the number of ANOVAs identified by both Y.J. and H.J., and the number of ANOVAs only identified by Y.J. or H.J. We are unable to evaluate this consistency as there is no information on the number of analyses that are not identified as ANOVA by both Y.J. and H.J. Therefore, we will only display the above numbers without calculating Cohen’s Kappa. Next, both Y.J. and H.J. re-reviewed these articles and the intra-rater consistency will be performed on the information after matching the numbers of identified ANOVAs for each article.

1.1.2.1 With mismatched numbers of ANOVAs

Load data for analysis information of the first 20 articles in each journal:

# read the coding sheet for the first round by YJ (before)
df_conv_YJ1_before <- 
  readxl::read_excel(here::here("data", "CodingSheet_first20.xlsx"),
                     sheet="AnalysisInfo_first20_YJ_before",
                     range="A1:AR205") %>%
  filter(!is.na(`Paper ID`)) %>% 
  select(`Paper ID`, First20, Exp_ID, DV) %>% 
  mutate(rater = "YJ")

# read the coding sheet for the first round by HJ (before)
df_conv_HJ1_before <- 
  readxl::read_excel(here::here("data", "CodingSheet_first20.xlsx"),
                     sheet="AnalysisInfo_first20_HJ_before",
                     range="A1:AR205") %>%
  filter(!is.na(`Paper ID`)) %>% 
  select(`Paper ID`, First20, Exp_ID, DV) %>% 
  mutate(rater = "HJ")

Contrast the number of ANOVAs identified by Y.J. and H.J.:

df_conv_before <- bind_rows(df_conv_YJ1_before, df_conv_HJ1_before) %>% 
  group_by(rater, `Paper ID`) %>% 
  summarize(count = n(), .groups = "drop") %>% 
  pivot_wider(names_from = rater, values_from = count, values_fill = 0) %>% 
  mutate(diff = HJ-YJ,
         YJ1_HJ1 = pmin(HJ, YJ),
         YJ0_HJ1 = ifelse(diff>0, diff, 0),
         YJ1_HJ0 = ifelse(diff<0, -diff, 0)) %>% 
  arrange(as.numeric(str_remove(`Paper ID`, "P"))) 

sum_conv_before <- df_conv_before %>% 
  summarize(YJ1_HJ1 = sum(YJ1_HJ1),
            YJ0_HJ1 = sum(YJ0_HJ1),
            YJ1_HJ0 = sum(YJ1_HJ0))
sum_conv_before 

Since it may happen that when one rater identified 10 ANOVAs and the other identified 8 ANOVAs, only 7 of them were identical. In other words, the first and the second raters missed 1 and 3 ANOVAs, respectively. We additionally check this case for all articles with different numbers of ANOVAs identified by Y.J. and H.J. The numbers of identified ANOVA after corrections are:

cmatrix_first20_before <- tibble(
  Rating = c("anova_HJ_1", "anova_HJ_0"),
  anova_YJ_1 = c(sum_conv_before$YJ1_HJ1-1, sum_conv_before$YJ1_HJ0+1),
  anova_YJ_0 = c(sum_conv_before$YJ0_HJ1+1, NA),
)
cmatrix_first20_before

1.1.2.2 After matching the number of ANOVAs in each article

Load analysis information after matching the number of ANOVAs in each article:

# read the coding sheet for the first round by YJ (after)
df_conv_YJ1_after <- 
  readxl::read_excel(here::here("data", "CodingSheet_first20.xlsx"),
                     sheet="AnalysisInfo_first20_YJ_after",
                     range="A1:AR205") %>%
  filter(!is.na(`Paper ID`)) %>% 
  select(-Title) %>% 
  mutate(row_idx = c(2:205))

# read the coding sheet for the first round by HJ (after)
df_conv_HJ1_after <- 
  readxl::read_excel(here::here("data", "CodingSheet_first20.xlsx"),
                     sheet="AnalysisInfo_first20_HJ_after",
                     range="A1:AR205") %>%
  filter(!is.na(`Paper ID`)) %>% 
  select(-Title) %>% 
  mutate(row_idx = c(2:205))

1.1.2.2.1 Experimental designs

Cohen’s Kappa for experimental designs:

# experimental designs
design_list <- c("N_factor", "max_No_levels")

kappa_ed <- lapply(design_list, 
                  thekappa, # available in R/func_kappa.R
                  df1=df_conv_YJ1_after, 
                  df2=df_conv_HJ1_after, 
                  isweighted=TRUE) %>%
  bind_rows()
kappa_ed

1.1.2.2.2 Convention A

Cohen’s Kappa for convention A:

# Convention A
con_a_list <- c("a1_report_all", "a3_correction")

kappa_a <- lapply(con_a_list, 
                  thekappa, # available in R/func_kappa.R
                  df1=filter(df_conv_YJ1_after, a2_multiway==1),
                  df2=filter(df_conv_HJ1_after, a2_multiway==1)) %>%
  bind_rows()
kappa_a

Confusion matrix for a1_report_all:

as.data.frame(kappa_a$conf_mat[[1]])

Rows with different information for a1_report_all:

bind_cols(a1_report_all_YJ = df_conv_YJ1_after$a1_report_all,
          df_conv_HJ1_after) %>%  
  select(a1_report_all_YJ, 
         a1_report_all_HJ = a1_report_all,
         row_idx, `Paper ID`, Exp_ID, DV, a2_multiway, everything()) %>% 
  filter(a1_report_all_YJ != a1_report_all_HJ,
         a2_multiway==1)

Confusion matrix for a3_correction:

as.data.frame(kappa_a$conf_mat[[2]])

1.1.2.2.3 Convention B

Main effect and post-hoc analysis.

Cohen’s Kappa for convention B:

# Convention B
con_b_wlist <- 
  c("b1_N_factor3+", # number of factors with three or more levels
    "b2_N_significant_b1", # number of (reported) significant factors with three or more levels
    "b3_N_non-sig_b1") # number of (reported) non-significant factors with three or more levels

kappa_bw <- lapply(con_b_wlist, 
                   thekappa, # available in R/func_kappa.R
                   isweighted=TRUE,
                   df1=filter(df_conv_YJ1_after, 
                              max_No_levels>2), 
                   df2=filter(df_conv_HJ1_after, 
                              max_No_levels>2)) %>%
  bind_rows()
kappa_bw

Confusion matrix for b2_N_significant_b1:

as.data.frame(kappa_bw$conf_mat[[2]])

Rows with different information on b2_N_significant_b1:

bind_cols(b2_N_significant_b1_YJ = df_conv_YJ1_after$`b2_N_significant_b1`,
          df_conv_HJ1_after) %>%  
  filter(max_No_levels>2) %>% 
  select(b2_N_significant_b1_YJ, 
         b2_N_significant_b1_HJ = `b2_N_significant_b1`,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(b2_N_significant_b1_YJ != b2_N_significant_b1_HJ)

Confusion matrix for b3_N_non-sig_b1:

as.data.frame(kappa_bw$conf_mat[[3]])

Rows with different information on b3_N_non-sig_b1:

bind_cols(`b3_N_non-sig_b1_YJ` = df_conv_YJ1_after$`b3_N_non-sig_b1`,
          df_conv_HJ1_after) %>%  
  filter(max_No_levels>2) %>% 
  select(`b3_N_non-sig_b1_YJ`, 
         `b3_N_non-sig_b1_HJ` = `b3_N_non-sig_b1`,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(`b3_N_non-sig_b1_YJ` != `b3_N_non-sig_b1_HJ`)

Only evaluate consistency of b7_main_name using the rows with consistent information in con_b_wlist:

tmp_exclude_b1 <- c(186, 187) # rows with inconsistent information
kappa_bu1 <- thekappa(
  df1=filter(df_conv_YJ1_after, 
             max_No_levels>2, 
             !(row_idx %in% tmp_exclude_b1)),
  df2=filter(df_conv_HJ1_after,
             max_No_levels>2,
             !(row_idx %in% tmp_exclude_b1)),
  colname = "b7_main_name") # which factor was checked (1 denotes the only one; NA denotes 'not applicable')

kappa_bu1

Confusion matrix for b7_main_name:

as.data.frame(kappa_bu1$conf_mat[[1]])

Rows with different information on b7_main_name:

bind_cols(`b7_main_name_YJ` = df_conv_YJ1_after$`b7_main_name`,
          df_conv_HJ1_after) %>%  
  filter(max_No_levels>2,
         !(row_idx %in% tmp_exclude_b1)) %>% 
  select(`b7_main_name_YJ`, 
         `b7_main_name_HJ` = `b7_main_name`,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(`b7_main_name_YJ` != `b7_main_name_HJ`)

Only evaluate consistency of b8_posthoc using the rows with consistent information in con_b_wlist and b7_main_name:

tmp_exclude_b2 <- c(186, 187, 199, 200) # rows with inconsistent information
kappa_bu2 <- thekappa(
  df1=filter(df_conv_YJ1_after, 
             max_No_levels>2, 
             b7_main_name != "NA",
             !(row_idx %in% tmp_exclude_b2)),
  df2=filter(df_conv_HJ1_after,
             max_No_levels>2,
             b7_main_name != "NA",
             !(row_idx %in% tmp_exclude_b2)),
  colname = "b8_posthoc") 

kappa_bu2

Confusion matrix for b8_posthoc:

as.data.frame(kappa_bu2$conf_mat[[1]])

Rows with different information on b8_posthoc:

bind_cols(b8_posthoc_YJ = df_conv_YJ1_after$b8_posthoc,
          df_conv_HJ1_after) %>%  
  filter(max_No_levels>2, 
         b7_main_name != "NA",
         !(row_idx %in% tmp_exclude_b2)) %>% 
  select(b8_posthoc_YJ, 
         b8_posthoc_HJ = b8_posthoc,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(b8_posthoc_YJ != b8_posthoc_HJ)

Only evaluate consistency of b9_correction using rows with consistent information in con_b_wlist and performed post-hoc analysis (b8_posthoc==1):

tmp_exclude_b3 <- c(174, 186, 187, 199, 200) # rows with inconsistent information
kappa_bu3 <- thekappa(
  df1=filter(df_conv_YJ1_after, 
             max_No_levels>2,
             !(row_idx %in% tmp_exclude_b3),
             b8_posthoc==1),
  df2=filter(df_conv_HJ1_after,
             max_No_levels>2,
             !(row_idx %in% tmp_exclude_b3),
             b8_posthoc==1),
  colname = "b9_correction") 

kappa_bu3

Confusion matrix for b9_correction:

as.data.frame(kappa_bu3$conf_mat[[1]])

Rows with different information on b9_correction:

bind_cols(b9_correction_YJ = df_conv_YJ1_after$b9_correction,
          df_conv_HJ1_after) %>%  
  filter(max_No_levels>2,
         !(row_idx %in% tmp_exclude_b3),
         b8_posthoc==1) %>% 
  select(b9_correction_YJ, 
         b9_correction_HJ = b9_correction,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(b9_correction_YJ != b9_correction_HJ)

Summary for Convention B:

(kappa_b <- bind_rows(kappa_bw, kappa_bu1, kappa_bu2, kappa_bu3))

1.1.2.2.4 Convention C

Interaction and simple effect analysis

Cohen’s Kappa for convention C:

# Convention C
con_c_wlist <- c("c1_N_2way_inter", "c2_N_sig_in_c1", "c3_N_non-sig_in_c1")

kappa_cw <- lapply(con_c_wlist, 
                   thekappa, # available in R/func_kappa.R
                   isweighted=TRUE,
                   df1=filter(df_conv_YJ1_after, N_factor>1),
                   df2=filter(df_conv_HJ1_after, N_factor>1)) %>%
  bind_rows()

kappa_cw

Confusion matrix for c1_N_2way_inter:

as.data.frame(kappa_cw$conf_mat[[1]])

Confusion matrix for c2_N_sig_in_c1:

as.data.frame(kappa_cw$conf_mat[[2]])

Rows with different information on c2_N_sig_in_c1:

bind_cols(c2_N_sig_in_c1_YJ = df_conv_YJ1_after$c2_N_sig_in_c1,
          df_conv_HJ1_after) %>%  
  filter(N_factor>1) %>% 
  select(c2_N_sig_in_c1_YJ, 
         c2_N_sig_in_c1_HJ = c2_N_sig_in_c1,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(c2_N_sig_in_c1_YJ != c2_N_sig_in_c1_HJ)

Confusion matrix for c3_N_non-sig_in_c1:

as.data.frame(kappa_cw$conf_mat[[3]])

Rows with different information on c3_N_non-sig_in_c1:

bind_cols(`c3_N_non-sig_in_c1_YJ` = df_conv_YJ1_after$`c3_N_non-sig_in_c1`,
          df_conv_HJ1_after) %>%  
  filter(N_factor>1) %>% 
  select(`c3_N_non-sig_in_c1_YJ`, 
         `c3_N_non-sig_in_c1_HJ` = `c3_N_non-sig_in_c1`,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(`c3_N_non-sig_in_c1_YJ` != `c3_N_non-sig_in_c1_HJ`)

Only evaluate consistency of c7_2way_name using the rows with consistent information in con_c_wlist:

kappa_cu1 <- thekappa(df1=filter(df_conv_YJ1_after, 
                                 N_factor>1),
                      df2=filter(df_conv_HJ1_after, 
                                 N_factor>1),
                      colname = "c7_2way_name")
kappa_cu1

Confusion matrix for c7_2way_name:

as.data.frame(kappa_cu1$conf_mat[[1]])

Rows with different information on c7_2way_name:

diff_c7 <- bind_cols(c7_2way_name_YJ = df_conv_YJ1_after$c7_2way_name,
                     df_conv_HJ1_after) %>%  
  filter(N_factor > 1) %>% 
  select(c7_2way_name_YJ, 
         c7_2way_name_HJ = c7_2way_name,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(tolower(c7_2way_name_YJ) != tolower(c7_2way_name_HJ))
diff_c7

Only evaluate consistency of c8_simple_effect using the rows with consistent information in con_c_wlist and checked main effects (c7_2way_name != "NA"):

tmp_exclude_c2 <- diff_c7$row_idx

kappa_cu2 <- thekappa(
  df1=filter(df_conv_YJ1_after, 
             N_factor>1,
             c7_2way_name != "NA",
             !(row_idx %in% tmp_exclude_c2)),
  df2=filter(df_conv_HJ1_after, 
             N_factor>1,
             c7_2way_name != "NA",
             !(row_idx %in% tmp_exclude_c2)),
  colname = "c8_simple_effect")

kappa_cu2

Confusion matrix for c8_simple_effect:

as.data.frame(kappa_cu2$conf_mat[[1]])

Rows with different information on c8_simple_effect:

diff_c8 <- bind_cols(c8_simple_effect_YJ = df_conv_YJ1_after$c8_simple_effect,
                     df_conv_HJ1_after) %>%  
  filter(N_factor>1,
         c7_2way_name != "NA",
         !(row_idx %in% tmp_exclude_c2)) %>% 
  select(c8_simple_effect_YJ, 
         c8_simple_effect_HJ = c8_simple_effect,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(c8_simple_effect_YJ != c8_simple_effect_HJ)

diff_c8

Only evaluate consistency of c9_correction using the rows with consistent information in con_c_wlist and performed simple effect analysis (c8_simple_effect==1):

tmp_exclude_c3 <- c(diff_c7$row_idx, diff_c8$row_idx)

kappa_cu3 <- thekappa(
  df1=filter(df_conv_YJ1_after,
             N_factor>1,
             c8_simple_effect==1,
             !(row_idx %in% tmp_exclude_c3)),
  df2=filter(df_conv_HJ1_after, 
             N_factor>1,
             c8_simple_effect==1,
             !(row_idx %in% tmp_exclude_c2)),
  colname = "c9_correction")

kappa_cu3

Confusion matrix for c9_correction:

as.data.frame(kappa_cu3$conf_mat[[1]])

Rows with different information on c9_correction:

bind_cols(c9_correction_YJ = df_conv_YJ1_after$c9_correction,
          df_conv_HJ1_after) %>%  
  filter(N_factor>1,
         c8_simple_effect==1,
         !(row_idx %in% tmp_exclude_c2)) %>% 
  select(c9_correction_YJ, 
         c9_correction_HJ = c9_correction,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(c9_correction_YJ != c9_correction_HJ) 

Summary for Convention C:

(kappa_c <- bind_rows(kappa_cw, kappa_cu1, kappa_cu2, kappa_cu3))

1.1.2.2.5 Convention D and E

Separate or collapse factors and perform more ANOVAs depending on the statistical significance of a higher-level interaction.

Cohen’s Kappa for the significance of the highest-level interaction:

# Convention D & E
kappa_d1e1 <- thekappa(df1=filter(df_conv_YJ1_after, 
                                  N_factor>2),
                       df2=filter(df_conv_HJ1_after, 
                                  N_factor>2),
                       colname="d1_e1_highest_sig")
kappa_d1e1

Confusion matrix for d1_e1_highest_sig:

as.data.frame(kappa_d1e1$conf_mat[[1]])

Rows with different information for d1_e1_highest_sig:

bind_cols(d1_e1_highest_sig_YJ = df_conv_YJ1_after$d1_e1_highest_sig,
          df_conv_HJ1_after) %>%  
  filter(N_factor > 2) %>% 
  select(d1_e1_highest_sig_YJ, 
         d1_e1_highest_sig_HJ = d1_e1_highest_sig,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(d1_e1_highest_sig_YJ != d1_e1_highest_sig_HJ)

Cohen’s Kappa for Convention D and E:

# separate one factor (for significant higher-level interaction)
kappa_d <- thekappa(df1=filter(df_conv_YJ1_after, 
                               N_factor>2, 
                               d1_e1_highest_sig==1,
                               !(row_idx %in% 135)),
                    df2=filter(df_conv_HJ1_after, 
                               N_factor>2, 
                               d1_e1_highest_sig==1,
                               !(row_idx %in% 135)),
                    colname="d1_separate")

# combine/collapse one factor (for non-significant higher-level interaction)
kappa_e <- thekappa(df1=filter(df_conv_YJ1_after, 
                               N_factor>2, 
                               d1_e1_highest_sig==0,
                               !(row_idx %in% 135)),
                    df2=filter(df_conv_HJ1_after, 
                               N_factor>2, 
                               d1_e1_highest_sig==0,
                               !(row_idx %in% 135)),
                    colname="e1_combine")

bind_rows(kappa_d, kappa_e)

Confusion matrix for d1_separate:

as.data.frame(kappa_d$conf_mat[[1]])

Confusion matrix for e1_combine:

as.data.frame(kappa_e$conf_mat[[1]])

Rows with different information for e1_combine:

bind_cols(e1_combine_YJ = df_conv_YJ1_after$e1_combine,
          df_conv_HJ1_after) %>%  
  filter(N_factor > 2,
         d1_e1_highest_sig==0,
         !(row_idx %in% 135)) %>% 
  select(e1_combine_YJ, 
         e1_combine_HJ = e1_combine,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(e1_combine_YJ != e1_combine_HJ)

Summary for Convnetion D and E:

(kappa_de <- bind_rows(kappa_d1e1, kappa_d, kappa_e))

1.1.2.2.6 Other information

Cohen’s Kappa for simple interaction analysis:

kappa_o1 <- thekappa(df1=df_conv_YJ1_after,
                     df2=df_conv_HJ1_after,
                     colname="Simple interaction analysis (appliable?)")

kappa_o2 <- thekappa(df1=filter(df_conv_YJ1_after, 
                                `Simple interaction analysis (appliable?)`==1),
                     df2=filter(df_conv_HJ1_after, 
                                `Simple interaction analysis (appliable?)`==1),
                     colname="Simple interaction analysis (applied?)")

bind_rows(kappa_o1, kappa_o2)

Confusion matrix for Simple interaction analysis (appliable?):

as.data.frame(kappa_o1$conf_mat[[1]])

Confusion matrix for Simple interaction analysis (applied?):

as.data.frame(kappa_o2$conf_mat[[1]])

Rows with different information for Simple interaction analysis (applied?) when it is applicable:

bind_cols(applied_YJ = df_conv_YJ1_after$`Simple interaction analysis (applied?)`,
          df_conv_HJ1_after) %>%  
  filter(`Simple interaction analysis (appliable?)`==1) %>% 
  select(applied_YJ, 
         applied_HJ = `Simple interaction analysis (applied?)`,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(applied_YJ != applied_HJ)

Cohen’s Kappa for pre-registration information:

olist <- c("Pre-registration", 
           "Registered reports")

kappa_o3 <- lapply(olist, 
                  thekappa, # available in R/func_kappa.R
                  df1=df_conv_YJ1_after,
                  df2=df_conv_HJ1_after) %>%
  bind_rows()

kappa_o3

Confusion matrix for Pre-registration:

as.data.frame(kappa_o3$conf_mat[[1]])

Rows with different information for Pre-registration:

bind_cols(`Pre-registration_YJ` = df_conv_YJ1_after$`Pre-registration`,
          df_conv_HJ1_after) %>%  
  select(`Pre-registration_YJ`, 
         `Pre-registration_HJ` = `Pre-registration`,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(`Pre-registration_YJ` != `Pre-registration_HJ`)

Confusion matrix for Registered reports:

as.data.frame(kappa_o3$conf_mat[[2]])

Rows with different information for Registered reports:

bind_cols(`Registered reports_YJ` = df_conv_YJ1_after$`Registered reports`,
          df_conv_HJ1_after) %>%  
  select(`Registered reports_YJ`, 
         `Registered reports_HJ` = `Registered reports`,
         row_idx, `Paper ID`, Exp_ID, DV, Design, everything()) %>% 
  filter(`Registered reports_YJ` != `Registered reports_HJ`)

Summary for other information:

(kappa_o <- bind_rows(kappa_o1, kappa_o2, kappa_o3))

1.1.2.2.7 Summary consistency for ANOVA conventions

(kappa_df <- bind_rows(kappa_ed, kappa_a, kappa_b, kappa_c, kappa_de, kappa_o))

Descriptive statistics of the consistency:

# the mean consistency
tibble(mean = mean(kappa_df$value),
       min = min(kappa_df$value),
       max = max(kappa_df$value))

1.2 Prevalence of ANOVA and the conventions

Load analysis information of all reviewed ANOVAs:

# read the coding sheet for the first round (after resolving inconsistency)
# df_conv_1 <- readxl::read_excel(here::here("data", "CodingSheet_first20.xlsx"),
#                                  sheet="AnalysisInfo_first20",
#                                  range="A1:AR205") %>%
#   select(-Title)
# # read the coding sheet for the second round by YJ
# df_conv_YJ2 <- readxl::read_excel(here::here("data", "CodingSheet_pilot.xlsx"), 
#                                  sheet="AnalysisInfo_2nd_YJ",
#                                  range="A1:AO39") %>% 
#   select(-Title) 
# # read the coding sheet for the second round by YJ
# df_conv_HJ2 <- readxl::read_excel(here::here("data", "CodingSheet_pilot.xlsx"), 
#                                  sheet="AnalysisInfo_2nd_YJ",
#                                  range="A1:AO39") %>% 
#   select(-Title) 
# 
# df_conv <- bind_rows(df_conv_1, df_conv_YJ2, df_conv_HJ2) %>% 
#   left_join(df_publ, by="Paper ID")

# df_conv <- readxl::read_excel(here::here("data", "CodingSheet_pilot.xlsx"), 
#                                  sheet="AnalysisInfo_HJ",
#                                  range="A1:AO39") %>% 
#   left_join(df_publ, by="Paper ID") %>% 
#   select(-Title) 

1.2.1 Screening information

Articles included at each stage:

# number of articles included at each stage
df_included_all <- df_mainscreen %>%  
  summarize(journal = "Overall",
            In1_all = n(),
            In2_title = sum(1-Title_Prescreening),
            In3_research = sum(`Research article`),
            In4_ANOVA = sum(`Use ANOVA`),
            In5_include = sum(Included))

# number of articles at different stages
df_included <- df_mainscreen %>% 
  group_by(journal) %>% 
  summarize(In1_all = n(),
            In2_title = sum(1-Title_Prescreening),
            In3_research = sum(`Research article`),
            In4_ANOVA = sum(`Use ANOVA`),
            In5_include = sum(Included)) %>% 
  add_row(df_included_all, .before=1)
df_included

The percentage of included articles:

# Percentage (relative to the total number of articles)
df_included %>% 
  transmute(journal, In1_all,
            In2_title = paste0(round((In2_title/In1_all)*100, 2), '%'),
            In3_research = paste0(round((In3_research/In1_all)*100, 2), '%'),
            In4_ANOVA = paste0(round((In4_ANOVA/In1_all)*100, 2), '%'),
            In5_include = paste0(round((In5_include/In1_all)*100, 2), '%'))

The number of excluded articles at each stage:

# number of articles excluded at different stages
df_excluded <- df_mainscreen %>% 
  group_by(journal) %>% 
  summarize(Ex1_all = n(),
            Ex2_title = sum(Title_Prescreening),
            Ex3_research = sum(1-`Research article`),
            Ex4_ANOVA = sum(1-`Use ANOVA`),
            Ex5_exclude = sum(1-Included)) %>% 
  add_row( # add summary row
    summarize(df_mainscreen, 
              journal = "Overall",
              Ex1_all = n(),
              Ex2_title = sum(Title_Prescreening),
              Ex3_research = sum(1-`Research article`),
              Ex4_ANOVA = sum(1-`Use ANOVA`),
              Ex5_exclude = sum(1-Included)),
    .before=1)
df_excluded

1.2.2 ANOVA prevalence

The prevalence of ANOVAs:

# The percentage of articles employing ANOVA among all research articles:  
df_included %>% 
  mutate(prevalence_anova = round(In4_ANOVA/In3_research*100, 2))

1.2.3 ANOVA complexity

# df_complex <- df_conv %>%
#   filter(included==1) %>%
#   transmute(`Paper ID`, N_factor, max_No_levels, Design,
#             Design_int = str_split(Design, "_"),
#             Max_level = sapply(Design_int, function(x){max(as.numeric(x))}),
#             Avg_No_levels = sapply(Design_int, function(x){mean(as.numeric(x))}),
#             Complexity = sapply(Design_int, function(x){prod(as.numeric(x))}),
#             journal)

Across ANOVAs:

# df_complex %>% 
#   group_by(journal) %>% 
#   summarize(mean_N_factor = mean(N_factor),
#             mean_max_N_levels = mean(max_No_levels),
#             mean_complexity = mean(Complexity)) %>% 
#   add_row( # add summary row
#     summarize(df_complex, 
#               journal = "All",
#               mean_N_factor = mean(N_factor),
#               mean_max_N_levels = mean(max_No_levels),
#               mean_complexity = mean(Complexity))) 

# plot_design <- df_complex %>% 
#   mutate(design_ = str_replace_all(Design, '_', '*')) %>% 
#   ggplot(aes(x = reorder(design_, Complexity))) +
#   geom_histogram(stat="count") +
#   labs(x = "Experimental designs")
# plot_design

alluvial figure (?) [left side is the experimental designs and the right side is the numbers of conditions.]

# df_complex %>% 
#   select(Design, Complexity) %>% 
#   group_by(Design, Complexity) %>% 
#   summarize(count = n(), .groups="drop") 

1.2.4 ANOVA convention prevalence

1.2.4.1 Convention A: inspect all effects

# df_conv_a <- df_conv %>% 
#   filter(a2_multiway==1) %>% 
#   select(`Paper ID`, starts_with("a"), journal)
# 
# # percentage of ANOVAs reporting/inspecting all effects among multi-way ANOVA
# df_conv_a_all <- summarize(df_conv_a, 
#                            journal = "All",
#                            conA = mean(a1_report_all))
# 
# df_conv_a %>% 
#   group_by(journal) %>% 
#   summarize(conA = mean(a1_report_all)) %>% 
#   add_row(df_conv_a_all) # add summary row

# df_conv_a_adjust <- df_conv_a %>% 
#   group_by(a3_correction) %>% 
#   summarize(journal = "All", count = n())
# 
# df_conv_a %>% 
#   group_by(journal, a3_correction) %>% 
#   summarize(count = n(), .groups="drop") %>% 
#   add_row(df_conv_a_adjust) 

1.2.4.2 Convention B: post-hoc

# df_conv_b <- df_conv %>% 
#   filter(max_No_levels>2) %>% 
#   select(`Paper ID`, max_No_levels, starts_with("b"), journal)
# 
# # percentage of (not) performing post-hoc analysis when the main effect is (not) significant
# df_conv_b_artile <- df_conv_b %>% 
#   mutate(isConB = b5_main_effect_sig == b7_posthoc) %>% 
#   group_by(`Paper ID`, journal) %>% 
#   summarize(meanB_art = mean(isConB), .groups="drop")
# 
# # show the result
# df_conv_b_all <- df_conv_b_artile %>% 
#   summarize(journal = "All", conB = mean(meanB_art))
# 
# df_conv_b_artile %>% 
#   group_by(journal) %>% 
#   summarize(conB = mean(meanB_art)) %>% 
#   add_row(df_conv_b_all) # add summary row

Multiple comparison corrections (per ANOVA)

# df_conv_b_adjust <- df_conv_b %>% 
#   filter(b5_main_effect_sig == 1) %>% 
#   group_by(b8_correction) %>% 
#   summarize(journal = 'All', count = n())
#   
# df_conv_b %>% 
#   filter(b5_main_effect_sig == 1) %>% 
#   group_by(journal, b8_correction) %>% 
#   summarize(count = n(), .groups="drop") %>% 
#   add_row(df_conv_b_adjust) 

1.2.4.3 Convention C: simple effect analysis

# df_conv_c <- df_conv %>% 
#   filter(N_factor>1) %>% 
#   select(`Paper ID`, N_factor, starts_with("c") & !starts_with("Con"), journal)
# 
# # percentage of (not) performing simple analysis when the interaction is (not) significant
# df_conv_c_artile <- df_conv_c %>% 
#   mutate(isConC = c5_inter_sig == c7_simple_effect) %>% 
#   group_by(`Paper ID`, journal) %>% 
#   summarize(meanC_art = mean(isConC), .groups="drop")
# 
# # show the result
# df_conv_c_all <- df_conv_c_artile %>% 
#   summarize(journal = "All", conC = mean(meanC_art))
# 
# df_conv_c_artile %>% 
#   group_by(journal) %>% 
#   summarize(conC = mean(meanC_art)) %>% 
#   add_row(df_conv_c_all) # add summary row

Multiple comparison corrections (per ANOVA)

# df_conv_c_adjust <- df_conv_c %>% 
#   filter(c5_inter_sig == 1) %>% 
#   group_by(c8_correction) %>% 
#   summarize(journal = 'All', count = n())
#   
# df_conv_c %>% 
#   filter(c5_inter_sig == 1) %>% 
#   group_by(journal, c8_correction) %>% 
#   summarize(count = n(), .groups="drop") %>% 
#   add_row(df_conv_c_adjust)

1.2.4.4 Convention D: separate

# df_conv_de <- df_conv %>% 
#   filter(N_factor>2) %>% 
#   select(`Paper ID`, N_factor, d1_e1_highest_sig, d1_separate, e1_combine, journal)
# 
# # 
# df_conv_d_artile <- df_conv_de %>% 
#   filter(d1_e1_highest_sig == 1) %>% 
#   group_by(`Paper ID`, journal) %>% 
#   summarize(meanD_art = mean(d1_separate), .groups="drop")
# 
# # show the result
# df_conv_d_all <- df_conv_d_artile %>% 
#   summarize(journal = "All", conD = mean(meanD_art))
# 
# df_conv_d_artile %>% 
#   group_by(journal) %>% 
#   summarize(conD = mean(meanD_art)) %>% 
#   add_row(df_conv_d_all) # add summary row

1.2.4.5 Convention E: combine

# # 
# df_conv_e_artile <- df_conv_de %>% 
#   filter(d1_e1_highest_sig == 0) %>% 
#   group_by(`Paper ID`, journal) %>% 
#   summarize(meanE_art = mean(e1_combine), .groups="drop")
# 
# # show the result
# df_conv_e_all <- df_conv_e_artile %>% 
#   summarize(journal = "All", conE = mean(meanE_art))
# 
# df_conv_e_artile %>% 
#   group_by(journal) %>% 
#   summarize(conE = mean(meanE_art)) %>% 
#   add_row(df_conv_e_all) # add summary row

2 Simulations for false positive rates in ANOVA conventions

Some texts.

2.1 Preparations

two_colors <- c("#D55E00", "#56B4E9")
custom_clr <- c("#7F7F7F", "#E64B35", "#56B4E9", "#DF8F44", "#008B45", "#925E9F", "#415384", "#000000")

2.2 Simulating Type I error in ANOVA

Nsim <- 10000 # number of simulation
alphas <- c(0.001, 0.005, 0.01, 0.05, (1:4)/10)  # 0.001 ~ 0.4

2.2.1 ANOVA omnibus F-test

p_omni <- sim_omnibus(N_subj = 30, iter = Nsim, n_core=16, 
                      file_cache = here::here("simulation", "p_omni.rds"),
                      N_IV = 2:5)

2.2.1.1 FWER and CER

# Type I error rates in different situations
df_TypeI_omni <- sig_omnibus(p_omni, 0.05) %>% 
  select(-c(alpha, effnames, p.value, sig, omniF)) %>% 
  distinct() %>% 
  group_by(N_IV, adjust) %>% 
  summarize(TypeI_omniF = mean(sig_omniF), 
            TypeI_any = mean(sig_any),
            `TypeI_omniF&any` = mean(`sig_omni&any`),
            .groups = "drop") 

df_TypeI_omni

# Conclusion-based Type I error rate (CBER)
df_omni_cber <- df_TypeI_omni %>% 
  pivot_longer(starts_with("TypeI_"), names_to = "effnames", values_to = "TypeI") %>% 
  mutate(adjust = if_else(effnames == "TypeI_omniF", "uncorrected", adjust),
         adjust = fct_relevel(adjust, "uncorrected"),
         effnames = factor(effnames, 
                           levels = c("TypeI_omniF", "TypeI_any", "TypeI_omniF&any")),
         effnames = fct_recode(effnames, `omnibus\nF-test` = "TypeI_omniF", 
                              `Inspect all\n(FWER)`= "TypeI_any", 
                              `omnibus & all\n(CER)`="TypeI_omniF&any")) %>% 
  distinct() 

fig_simu_a <- df_omni_cber %>% 
  filter(N_IV == 2) %>% # only show results for ANOVAs with two factors
  ggplot(aes(x=as.numeric(adjust), y = TypeI*100, fill=adjust)) + 
  geom_col(width=1) +   
  facet_grid(cols=vars(effnames), scales="free_x", space="free_x", switch="x")+
  ggh4x::facetted_pos_scales(
    list(
      scale_x_continuous(breaks = 1, expand = c(0, 0.65)),
      scale_x_continuous(breaks = 1.5, expand = c(0, 0.65)),
      scale_x_continuous(breaks = 1.5, expand = c(0, 0.65))
    )) + 
  scale_fill_manual(values = custom_clr[c(1,5)],
                    breaks = levels(df_omni_cber$adjust)) + #,  
  geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
  scale_y_continuous(limits = c(0, 24), expand = c(0, 0)) +
  labs(x = "ANOVAs with two factors", y = "False Positive Rate (%)", fill = "") +
  theme(legend.position = c(0.75, 0.6),
        axis.ticks.length = unit(.15, "cm"),
        axis.text.x = element_blank(),
        panel.spacing.x = unit(0, "pt"),
        strip.placement = "outside",
        axis.title.x = element_text(margin = margin(t = 0)))
  
# ggsave(file.path("figures","fig_simu_TypeI_omni.png"), fig_simu_a, width = 4, height = 5.6)
fig_simu_a

# show results for ANOVAs with different factor numbers
for (i in unique(df_omni_cber$N_IV)) {
  
  nam <- paste("fig_simu_a_tmp", i, sep = "")
  
  assign(
    nam, df_omni_cber %>% 
      filter(N_IV == i) %>% 
      ggplot(aes(x=as.numeric(adjust), y = TypeI*100, fill=adjust)) + 
      geom_col(width=1) +   
      facet_grid(cols=vars(effnames), scales="free_x", space="free_x", switch="x")+
      ggh4x::facetted_pos_scales(
        list(
          scale_x_continuous(breaks = 1, expand = c(0, 0.65)),
          scale_x_continuous(breaks = 1.5, expand = c(0, 0.65)),
          scale_x_continuous(breaks = 1.5, expand = c(0, 0.65))
        )) + 
      scale_fill_manual(values = custom_clr[c(1,5)],
                        breaks = levels(df_omni_cber$adjust)) + #,  
      geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
      scale_y_continuous(limits = c(0, 84), expand = c(0, 0)) +
      labs(x = paste0("ANOVAs with ", i, " factors"), y = "False Positive Rate (%)", fill = "") +
      theme(legend.position = c(0.76, 0.6),
            axis.ticks.length = unit(.15, "cm"),
            axis.text.x = element_blank(),
            panel.spacing.x = unit(0, "pt"),
            strip.placement = "outside",
            axis.title.x = element_text(margin = margin(t = 0)))
  )
}

fig_simu_a_supp <- ggarrange(fig_simu_a_tmp2, fig_simu_a_tmp3, 
                             fig_simu_a_tmp4, fig_simu_a_tmp5, 
                             nrow = 1)
  
# ggsave(file.path("figures","fig_simu_TypeI_omni_supp.png"), fig_simu_a_supp, width = 14, height = 5.6)
fig_simu_a_supp

2.2.1.2 Single tests

# with omnibus F-test and significant effect (uncorrected)
df_sig_single1 <- sig_omnibus(p_omni, 0.05)  %>% 
  filter(adjust == "uncorrected") %>% 
  transmute(Method = "omnibus F-test",
            N_IV, effnames, 
            sig_cber = `sig_omni&any`,
            sig_cber_single = sig_cber & sig)

df_TypeI_single <- sig_omnibus(p_omni, 0.05)  %>% 
  filter(adjust == "Bonferroni") %>% 
  transmute(Method = "Bonferroni",
            N_IV, effnames,
            sig_cber = sig_any,
            sig_cber_single = sig_cber & sig) %>% 
  # combine the two df and make it to long format
  bind_rows(df_sig_single1) %>% 
  group_by(N_IV, effnames, Method) %>% 
  summarize(TypeI_cber = mean(sig_cber), 
            TypeI_cber_single = mean(sig_cber_single),
            .groups = "drop") %>% 
  pivot_wider(c(N_IV, Method, TypeI_cber), 
              names_from = effnames, values_from = TypeI_cber_single) %>% 
  pivot_longer(!c("N_IV", "Method"), 
               names_to = "effnames", values_to = "TypeI") %>% 
  filter(!is.na(TypeI)) %>% 
  # update the level names and order
  mutate(Method = factor(Method, levels=c("omnibus F-test", "Bonferroni")))

# Type I error rate for single tests (when the CBER is controlled)
df_TypeI_single %>%
  pivot_wider(names_from = effnames, values_from = TypeI) %>% 
  select(Method, N_IV, TypeI_cber, everything())

fig_simu_single_a <- df_TypeI_single %>% 
  filter(N_IV == 2) %>% 
  select(-c(N_IV)) %>% 
  mutate(effnames = factor(effnames),
         effnames = fct_relevel(effnames, "TypeI_cber", "A", "B"),
         effnames = fct_recode(effnames, CBER="TypeI_cber", `main effect A`="A",
                               `main effect B`="B", interaction="A:B"),
         Method = paste0(Method, "\n"),
         Method = factor(Method, levels=c("omnibus F-test\n", "Bonferroni\n"))) %>% 
  ggplot(aes(x = Method, y = TypeI*100, fill = effnames, group = effnames)) +
  geom_col(position = "dodge") +
  # facet_grid(~ omni, switch = "x") + 
  scale_fill_manual(values = custom_clr[1:4]) +
  geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
  scale_y_continuous(limits = c(0, 7), expand = c(0, 0), breaks = 1:6) +
  labs(x = "ANOVAs with two factors", y = "False Positive Rate (%)", fill = "") +
  guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
  theme(legend.position = c(0.55, 0.87),
        legend.direction="horizontal",
        axis.ticks.length = unit(.15, "cm"),
        strip.placement = "outside",
        axis.title.x = element_text(margin = margin(t = 5)))

# ggsave(file.path("figures","fig_simu_single_omni.png"), fig_simu_single_a, width = 4, height = 5.6)
fig_simu_single_a

df_TypeI_single_plot <- df_TypeI_single %>% 
  mutate(effnames = factor(effnames),
         effnames = fct_relevel(effnames, "TypeI_cber"),
         effnames = fct_recode(effnames, CBER="TypeI_cber"),
         Method = as_factor(Method)) 

for (i in unique(df_TypeI_single_plot$N_IV)) {
  
  nam <- paste("fig_simu_single_a_supp_tmp", i, sep = "")
  
  assign(
    nam, df_TypeI_single_plot %>% 
      filter(N_IV==i) %>% 
      ggplot(aes(x = Method, y = TypeI*100, fill = effnames, group = effnames)) +
      geom_col(position = "dodge") +
      # facet_wrap(vars(N_IV), switch = "x") +
      scale_fill_manual(values = c(custom_clr[c(1,3)], brewer.pal(8,"Set1"), 
                                   brewer.pal(7, "Set2"), brewer.pal(8,"Set3"), brewer.pal(7,"Accent"))) + 
      geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
      scale_y_continuous(limits = c(0, 7), expand = c(0, 0), breaks = 1:6) +
      labs(x = paste0("ANOVAs with ", as.character(i), " factors"), y = "False Positive Rate (%)", fill = "") +
      guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
      theme(legend.position = "none",
            # legend.direction="horizontal",
            axis.ticks.length = unit(.15, "cm"),
            strip.placement = "outside",
            axis.title.x = element_text(margin = margin(t = 5)))
  )
}

fig_simu_single_a_supp <- ggarrange(
  ggarrange(fig_simu_single_a_supp_tmp2, fig_simu_single_a_supp_tmp5,
            ncol=2, widths = c(0.3, 0.7)), 
  ggarrange(fig_simu_single_a_supp_tmp3, fig_simu_single_a_supp_tmp4,
            ncol=2, widths = c(0.4, 0.6)),
  nrow = 2)

# ggsave(file.path("figures","fig_simu_single_omni_supp.png"),
#        fig_simu_single_a_supp, width = 10, height = 7)
fig_simu_single_a_supp

2.2.2 Main effect and post-hoc analysis

# sim_main_posthoc is available in R/funcs.R
p_main_posthoc <- sim_main_posthoc(N_subj = 30, iter = Nsim, n_core=16, 
                                   file_cache = here::here("simulation", "p_main_posthoc.rds"),
                                   N_levels = c(3,4,5,6)) 

2.2.2.1 FWER and CER

df_TypeI_posthoc <- sig_main_posthoc(p_main_posthoc, 0.05) %>% 
  group_by(N_level, alpha, adjust) %>% 
  summarize(TypeI_main = mean(sig_main),
            TypeI_postonly = mean(sig_anypost),
            `TypeI_main&posthoc` = mean(`sig_main&posthoc`), 
            .groups = "drop")

df_TypeI_posthoc

Claim significant results only when both the main effect and at least one of the post-hoc tests (with different multiple comparison corrections) are significant:

df_TypeI_posthoc %>% 
  filter(adjust != "uncorrected") %>% # only when using multiple comparison correction
  select(-c(TypeI_main, TypeI_postonly)) %>% 
  pivot_wider(c(N_level, alpha), 
              names_from=adjust, values_from=`TypeI_main&posthoc`)

Claim significant results only when the main effect and at least one of the post-hoc tests (without multiple comparison corrections) are significant:

df_TypeI_posthoc %>% 
  filter(adjust == "uncorrected") %>%  # only when NOT using multiple comparison correction
  select(-c(TypeI_main, TypeI_postonly)) %>% 
  pivot_wider(c(N_level, alpha), 
              names_from=adjust, values_from=`TypeI_main&posthoc`)

Claim significant results only when at least one of the post-hoc tests (with different multiple comparison corrections) are significant (regardless of the main effect results):

df_TypeI_posthoc %>% 
  filter(adjust != "uncorrected") %>% # only when using multiple comparison correction
  select(-c(TypeI_main, `TypeI_main&posthoc`)) %>% 
  pivot_wider(c(N_level, alpha), 
              names_from=adjust, values_from=TypeI_postonly)

df_posthoc_tmp <- df_TypeI_posthoc %>% 
  pivot_longer(starts_with("TypeI_"), names_to = "effects", values_to = "TypeI") %>% 
  mutate(adjust = if_else(effects == "TypeI_main", "uncorrected", 
                          if_else(adjust == "uncorrected", "uncorrected", 
                                  str_to_title(as.character(adjust)))),
         adjust = fct_reorder(adjust, TypeI, mean, .desc=TRUE),
         adjust = fct_relevel(adjust, "uncorrected"),
         effects = factor(effects, 
                          levels = c("TypeI_main", "TypeI_postonly", "TypeI_main&posthoc")),
         effects = fct_recode(effects, `Main effect\n` ="TypeI_main", 
                              `Pairwise comparisons\n(six tests, FWER)`="TypeI_postonly", 
                              `Main effect & Pairwise comparisons\nwith six tests (CER)`="TypeI_main&posthoc")) %>% 
  distinct() 

fig_simu_b <- df_posthoc_tmp %>% 
  filter(N_level == 4) %>% 
  ggplot(aes(x=as.numeric(adjust), y = TypeI*100, fill=adjust)) + 
  geom_col(width=1) +   
  facet_grid(cols=vars(effects), scales="free_x", space="free_x", switch="x")+
  ggh4x::facetted_pos_scales(
    list(
      scale_x_continuous(breaks = 1, expand = c(0, 0.65)),
      scale_x_continuous(breaks = 3.5, expand = c(0, 0.65)),
      scale_x_continuous(breaks = 3.5, expand = c(0, 0.65))
    )) + 
  scale_fill_manual(values = custom_clr,
                    breaks = levels(df_posthoc_tmp$adjust)) + #,  
  geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
  scale_y_continuous(limits = c(0, 24), expand = c(0, 0)) +
  labs(x = "One-way ANOVAs with four levels", y = "False Positive Rate (%)", fill = "") +
  theme(legend.position = c(0.74, 0.63),
        axis.ticks.length = unit(.15, "cm"),
        axis.text.x = element_blank(),
        panel.spacing.x = unit(0, "pt"),
        strip.placement = "outside",
        axis.title.x = element_text(margin = margin(t = 0)))
  
# ggsave(file.path("figures","fig_simu_TypeI_mainpost.png"), fig_simu_b, width = 8, height = 5.6)
fig_simu_b

for (i in unique(df_posthoc_tmp$N_level)) {
  
  nam <- paste("fig_simu_b_tmp", i, sep = "")
  
  assign(
    nam, df_posthoc_tmp %>% 
    filter(N_level == i) %>% 
    ggplot(aes(x=as.numeric(adjust), y = TypeI*100, fill=adjust)) + 
    geom_col(width=1) +   
    facet_grid(cols=vars(effects), scales="free_x", space="free_x", switch="x")+
    ggh4x::facetted_pos_scales(
      list(
        scale_x_continuous(breaks = 1, expand = c(0, 0.65)),
        scale_x_continuous(breaks = 3.5, expand = c(0, 0.65)),
        scale_x_continuous(breaks = 3.5, expand = c(0, 0.65))
      )) + 
    scale_fill_manual(values = custom_clr,
                      breaks = levels(df_posthoc_tmp$adjust)) + #,  
    geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
    scale_y_continuous(limits = c(0, 45), expand = c(0, 0)) +
    labs(x = paste0("One-way ANOVAs with ", i, " levels"), y = "False Positive Rate (%)", fill = "") +
    theme(legend.position = c(0.75, 0.6),
          axis.ticks.length = unit(.15, "cm"),
          axis.text.x = element_blank(),
          panel.spacing.x = unit(0, "pt"),
          strip.placement = "outside",
          axis.title.x = element_text(margin = margin(t = 0)))
  )
  
}

fig_simu_b_supp <- ggarrange(fig_simu_b_tmp3, fig_simu_b_tmp4, 
                             fig_simu_b_tmp5, fig_simu_b_tmp6, 
                             nrow = 2, ncol = 2)

# ggsave(file.path("figures","fig_simu_TypeI_mainpost_supp.png"), fig_simu_b_supp, width = 13, height = 8)
fig_simu_b_supp

2.2.2.2 Single tests

# with omnibus F-test and significant effect
df_sig_main_single1 <- sig_main_posthoc(p_main_posthoc, 0.05) %>% 
  filter(adjust == "uncorrected") %>% 
  transmute(Method = "Main effect & Pairwise",
            N_level, contrast,
            sig_cber = `sig_main&posthoc`, # main and post-hoc (uncorrected)
            sig_cber_single = sig_cber & sig)

df_TypeI_main_single <- sig_main_posthoc(p_main_posthoc, 0.05) %>% 
  filter(adjust != "uncorrected") %>% 
  transmute(Method = str_to_title(adjust),
            N_level, contrast,
            sig_cber = `sig_anypost`, # post-hoc only (corrected)
            sig_cber_single = sig_cber & sig) %>% 
  # combine two df and make it to long format
  bind_rows(df_sig_main_single1) %>% 
  group_by(N_level, contrast, Method) %>% 
  summarize(TypeI_cber = mean(sig_cber), # correction and at least one sig
            TypeI_cber_single = mean(sig_cber_single),
            .groups = "drop") %>% 
  pivot_wider(c(N_level, Method, TypeI_cber), 
              names_from = contrast, values_from = TypeI_cber_single) %>% 
  pivot_longer(!c("N_level", "Method"), 
               names_to = "contrast", values_to = "TypeI") %>% 
  filter(!is.na(TypeI))
  
# Type I error rate for single tests (when the FWER is controlled)
df_TypeI_main_single %>%
  pivot_wider(names_from = contrast, values_from = TypeI) %>% 
  select(Method, N_level, TypeI_cber, everything())

fig_simu_single_b <- df_TypeI_main_single %>% 
  filter(N_level == 4) %>% 
  select(-c(N_level)) %>% 
  mutate(contrast = factor(contrast),
         contrast = fct_relevel(contrast, "TypeI_cber"),
         contrast = fct_recode(contrast, CBER="TypeI_cber"),
         contrast = str_replace_all(contrast, "IV", "level"),
         Method = if_else(str_starts(Method, "Main"),
                          paste0(Method, "\n(uncorrected)"),
                          paste0(Method, "\n(pairwise only)")),
         Method = as_factor(Method),
         Method = fct_reorder(Method, TypeI, mean, .desc=TRUE)) %>%
  ggplot(aes(x = Method, y = TypeI*100, fill = contrast, group = contrast)) +
  geom_col(position = "dodge") +
  scale_fill_manual(values = custom_clr) +
  geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
  scale_y_continuous(limits = c(0, 7), expand = c(0, 0), breaks = 1:6) +
  labs(x = "One-way ANOVAs with four levels", y = "False Positive Rate (%)", fill = "") +
  guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
  theme(legend.position = c(0.5, 0.87),
        legend.direction="horizontal",
        axis.ticks.length = unit(.15, "cm"),
        strip.placement = "outside",
        axis.title.x = element_text(margin = margin(t = 5)))

# ggsave(file.path("figures","fig_simu_single_mainpost.png"), 
#        fig_simu_single_b, width = 12, height = 4.8)
fig_simu_single_b

for (i in unique(df_TypeI_main_single$N_level)) {
  
  nam <- paste("fig_simu_single_b_tmp", i, sep = "")
  
  assign(
    nam, df_TypeI_main_single %>% 
      filter(N_level == i) %>% 
      select(-c(N_level)) %>% 
      mutate(contrast = factor(contrast),
             contrast = fct_relevel(contrast, "TypeI_cber"),
             contrast = fct_recode(contrast, CBER="TypeI_cber"),
             contrast = str_replace_all(contrast, "IV", "level"),
             Method = if_else(str_starts(Method, "Main"), 
                          paste0(Method, "\n(uncorrected)"),
                          paste0(Method, "\n(post-hoc only)")),
             Method = as_factor(Method),
             Method = fct_reorder(Method, TypeI, mean, .desc=TRUE)) %>%
      ggplot(aes(x = Method, y = TypeI*100, fill = contrast, group = contrast)) +
      geom_col(position = "dodge") +
      scale_fill_manual(values = c(custom_clr[c(1,3)], brewer.pal(8,"Set1"),
                                   brewer.pal(6,"Set3"))) +
      geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
      scale_y_continuous(limits = c(0, 7), expand = c(0, 0), breaks = 1:6) +
      labs(x = paste0("One-way ANOVAs with ", i, " levels"), y = "False Positive Rate (%)", fill = "") +
      guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
      theme(legend.position = c(0.5, 0.88),
            legend.direction="horizontal",
            axis.ticks.length = unit(.15, "cm"),
            strip.placement = "outside")
  )
}

fig_simu_single_b_supp <- ggarrange(fig_simu_single_b_tmp3, 
                                    fig_simu_single_b_tmp4, 
                                    ncol=2) %>% 
  ggarrange(fig_simu_single_b_tmp5,
            fig_simu_single_b_tmp6,
            ncol=1)

# ggsave(file.path("figures","fig_simu_single_mainpost_supp.png"), 
#        fig_simu_single_b_supp, width = 14, height = 14)

fig_simu_single_b_supp

2.2.3 Interaction and simple effect analysis

p_inter_simple <- sim_inter_simple(N_subj = 30, iter = Nsim, n_core=16, 
                                   file_cache = here::here("simulation", "p_inter_simple.rds"))

2.2.3.1 FWER and CER

Type I error rates for main effects and interaction (when different correction methods were applied to the simple effect analysis):

p_inter_simple %>% 
  select(iter, adjust, p_mainA, p_mainB, p_inter) %>% 
  distinct() %>% 
  group_by(adjust) %>% 
  summarize(mean(p_mainA<0.05), mean(p_mainB<0.05), mean(p_inter<0.05)) 

Type I error rates (for different combinations of simple effects and interaction):

df_TypeI_inter <- sig_inter_simple(p_inter_simple, 0.05) %>% 
  select(-c(contrast, sig_simple)) %>% 
  distinct() %>% 
  group_by(adjust, group) %>% 
  summarize(TypeI_inter = mean(sig_inter),
            TypeI_any_simple = mean(sig_any_simple),
            `TypeI_inter&any` = mean(`sig_inter&any`),
            .groups="drop") %>% 
  pivot_longer(c(TypeI_any_simple, `TypeI_inter&any`), 
               names_to = "effects", values_to = "TypeI") %>% 
  mutate(effects = paste(effects, group, sep="_")) 

# Type I error rates
df_TypeI_inter %>% 
  select(-group) %>% 
  pivot_wider(names_from = effects, values_from = TypeI) %>% 
  select(adjust, TypeI_inter, TypeI_any_simple_byAB_4, 
         `TypeI_inter&any_byAB_4`, everything()) %>% 
  filter(adjust != "uncorrected") 

df_TypeI_inter %>% 
  select(-group) %>% 
  filter(adjust == "uncorrected") 

df_inter_simp_tmp <- df_TypeI_inter %>% 
  filter(group != "byB_2") %>% 
  select(-c(group, TypeI_inter)) %>% 
  add_row(adjust="uncorrected", effects = "Interaction\n", 
          TypeI = unique(df_TypeI_inter$TypeI_inter)) %>% 
  mutate(
    adjust = if_else(adjust == "uncorrected", "uncorrected",
                     str_to_title(as.character(adjust))),
    adjust = fct_reorder(adjust, TypeI, mean, .desc=TRUE),
    adjust = fct_relevel(adjust, "uncorrected"),
    effects = factor(effects, 
                     levels = c("Interaction\n", 
                                "TypeI_any_simple_byA_2", "TypeI_any_simple_byAB_4",
                                "TypeI_inter&any_byA_2", "TypeI_inter&any_byAB_4")),
    effects = fct_recode(
      effects, 
      `Simple effect analysis\n(two tests, FWER)`= "TypeI_any_simple_byA_2", 
      `Simple effect analysis\n(four tests, FWER)`="TypeI_any_simple_byAB_4",
      `Interaction & simple effects\nwith two tests (CER)`="TypeI_inter&any_byA_2",
      `Interaction & simple effects\nwith four tests (CER)`= "TypeI_inter&any_byAB_4"))

fig_simu_c <- df_inter_simp_tmp %>% 
  ggplot(aes(x=as.numeric(adjust), y = TypeI*100, fill=adjust)) + 
  geom_col(width=1) +   
  facet_grid(cols=vars(effects), scales="free_x", space="free_x", switch="x") +
  ggh4x::facetted_pos_scales(
    list(
      scale_x_continuous(breaks = 1, expand = c(0, 0.65)),
      scale_x_continuous(breaks = 3, expand = c(0, 0.65)),
      scale_x_continuous(breaks = 3, expand = c(0, 0.65)),
      scale_x_continuous(breaks = 3, expand = c(0, 0.65)),
      scale_x_continuous(breaks = 3, expand = c(0, 0.65))
    )) + 
  scale_fill_manual(values = custom_clr[c(1,3:6)],
                    breaks = levels(df_inter_simp_tmp$adjust)) + #,
  geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
  scale_y_continuous(limits = c(0, 18), expand = c(0, 0)) +
  labs(x = "ANOVAs with 2 × 2 designs", y = "False Positive Rate (%)", fill = "") +
  theme(legend.position = c(0.84, 0.65),
        axis.ticks.length = unit(.15, "cm"),
        axis.text.x = element_blank(),
        panel.spacing.x = unit(0, "pt"),
        strip.placement = "outside",
        axis.title.x = element_text(margin = margin(t = 0)))
  
# ggsave(file.path("figures","fig_simu_TypeI_inter.png"), fig_simu_c, width = 10, height = 4)
fig_simu_c

df_inter_simp_tmp_supp <- df_TypeI_inter %>% 
  filter(group != "byA_2") %>% # the only difference from previous figure
  select(-c(group, TypeI_inter)) %>% 
  add_row(adjust="uncorrected", effects = "Interaction\n", 
          TypeI = unique(df_TypeI_inter$TypeI_inter)) %>% 
  mutate(
    adjust = if_else(adjust == "uncorrected", "uncorrected",
                     str_to_title(as.character(adjust))),
    adjust = fct_reorder(adjust, TypeI, mean, .desc=TRUE),
    adjust = fct_relevel(adjust, "uncorrected"),
    effects = factor(effects, 
                     levels = c("Interaction\n", 
                                "TypeI_any_simple_byB_2", "TypeI_any_simple_byAB_4",
                                "TypeI_inter&any_byB_2", "TypeI_inter&any_byAB_4")),
    effects = fct_recode(
      effects, 
      `Simple effect analysis\n(two tests, FWER)`= "TypeI_any_simple_byB_2", 
      `Simple effect analysis\n(four tests, FWER)`="TypeI_any_simple_byAB_4",
      `Interaction & simple effects\nwith two tests (CER)`="TypeI_inter&any_byB_2",
      `Interaction & simple effects\nwith four tests (CER)`= "TypeI_inter&any_byAB_4"))

fig_simu_c_supp <- df_inter_simp_tmp_supp %>% 
  ggplot(aes(x=as.numeric(adjust), y = TypeI*100, fill=adjust)) + 
  geom_col(width=1) +   
  facet_grid(cols=vars(effects), scales="free_x", space="free_x", switch="x") +
  ggh4x::facetted_pos_scales(
    list(
      scale_x_continuous(breaks = 1, expand = c(0, 0.65)),
      scale_x_continuous(breaks = 3, expand = c(0, 0.65)),
      scale_x_continuous(breaks = 3, expand = c(0, 0.65)),
      scale_x_continuous(breaks = 3, expand = c(0, 0.65)),
      scale_x_continuous(breaks = 3, expand = c(0, 0.65))
    )) + 
  scale_fill_manual(values = custom_clr[c(1,3:6)],
                    breaks = levels(df_inter_simp_tmp$adjust)) + #,
  geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
  scale_y_continuous(limits = c(0, 18), expand = c(0, 0)) +
  labs(x = "ANOVAs with 2 × 2 designs", y = "False Positive Rate (%)", fill = "") +
  theme(legend.position = c(0.84, 0.65),
        axis.ticks.length = unit(.15, "cm"),
        axis.text.x = element_blank(),
        panel.spacing.x = unit(0, "pt"),
        strip.placement = "outside",
        axis.title.x = element_text(margin = margin(t = 0)))
  
# ggsave(file.path("figures","fig_simu_TypeI_inter_supp.png"), fig_simu_c_supp, width = 10, height = 4)
fig_simu_c_supp

2.2.3.2 Single tests

# with significant interaction and uncorrected simple effects
df_sig_inter_single1 <- sig_inter_simple(p_inter_simple, 0.05) %>% 
  filter(adjust == "uncorrected") %>% 
  transmute(Method = "Inter. & simple\n(uncorrected)",
            group, contrast, 
            sig_cber = `sig_inter&any`,
            sig_cber_single = sig_cber & sig_simple,
            sig_cber_inter = sig_cber & sig_inter) 

# with corrected simple effects
df_TypeI_inter_single <- sig_inter_simple(p_inter_simple, 0.05) %>% 
  filter(adjust != "uncorrected") %>% 
  transmute(Method = paste0(str_to_title(adjust), "\n(pairwise only)"),
            group, contrast, 
            sig_cber = sig_any_simple,
            sig_cber_single = sig_cber & sig_simple,
            sig_cber_inter = sig_cber & sig_inter) %>% 
  # combine the two df and make it to long format
  bind_rows(df_sig_inter_single1) %>% 
  group_by(group, contrast, Method) %>% 
  summarize(TypeI_cber = mean(sig_cber),
            TypeI_cber_single = mean(sig_cber_single),
            TypeI_cber_inter = mean(sig_cber_inter),
            .groups="drop") %>% 
  pivot_wider(c(group, Method, TypeI_cber, TypeI_cber_inter), 
              names_from = contrast, values_from = TypeI_cber_single) %>% 
  pivot_longer(!c("group", "Method"), 
               names_to = "contrast", values_to = "TypeI") %>% 
  filter(!is.na(TypeI)) %>% 
  # update the level names and order
  mutate(contrast = case_when(
    contrast == "a1 b1 - a1 b2" ~ "a1_(b1 - b2)",
    contrast == "a1 b1 - a2 b1" ~ "b1_(a1 - a2)",
    contrast == "a1 b2 - a2 b2" ~ "b2_(a1 - a2)",
    contrast == "a2 b1 - a2 b2" ~ "a2_(b1 - b2)",
    contrast == "TypeI_cber" ~ "CBER",
    contrast == "TypeI_cber_inter" ~ "Interaction",
    TRUE ~ contrast),
    contrast = as_factor(contrast),
    contrast = fct_relevel(contrast, "CBER", "Interaction",
                           "a1_(b1 - b2)", "a2_(b1 - b2)"),
    Method = fct_reorder(Method, TypeI, mean, .desc=TRUE),
    Method = fct_relevel(Method, "Inter. & simple\n(uncorrected)"),
    grouplabel = if_else(group=="byAB_4", 
                         "Simple effect analysis with four tests",
                         "Simple effect analysis with two tests"),
    grouplabel = as_factor(grouplabel),
    grouplabel = fct_relevel(grouplabel, "Simple effect analysis with two tests")
  ) %>% 
  filter(contrast!="Interaction")

df_TypeI_inter_single

fig_simu_single_c <- df_TypeI_inter_single %>% 
  filter(group != "byB_2") %>% 
  ggplot(aes(x = Method, y = TypeI*100, fill = contrast, group = contrast)) +
  geom_col(position = "dodge") +
  facet_grid(cols=vars(grouplabel), switch="x") +
  scale_fill_manual(values = custom_clr) +
  geom_hline(yintercept = 5, linetype= "dashed", color="gray25") +
  scale_y_continuous(limits = c(0, 7), expand = c(0, 0), breaks = 1:6) +
  labs(x = "ANOVAs with 2 × 2 designs", y = "False Positive Rate (%)", fill = "") +
  guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
  theme(legend.position = c(0.5, 0.87),
        legend.direction="horizontal",
        axis.ticks.length = unit(.15, "cm"),
        strip.placement = "outside",
        axis.title.x = element_text(margin = margin(t = 0)))

# ggsave(file.path("figures","fig_simu_single_inter.png"), fig_simu_single_c, width = 12, height = 5)
fig_simu_single_c

2.3 Publication

2.3.1 Figure 2

fig_simu <- ggarrange(fig_simu_a, fig_simu_b,
          nrow = 1, labels = c("a", "b"),
          widths = c(0.35, 0.65)) %>% 
  ggarrange(fig_simu_c,
            nrow = 2, labels = c("", "c"))

# ggsave(file.path("figures","fig_simu.png"), fig_simu, width = 12, height = 8.7)
fig_simu

fig_simu_single <- ggarrange(fig_simu_single_a, fig_simu_single_b,
          nrow = 1, labels = c("a", "b"),
          widths = c(0.3, 0.7)) %>% 
   ggarrange(fig_simu_single_c,
            nrow = 2, labels = c("", "c"))

# ggsave(file.path("figures","fig_simu_single.png"), fig_simu_single, width = 12, height = 10.3)
fig_simu_single

2.3.2 Table 1

2.4 Supplemenatry

fig_simu_a_supp

ggsave(file.path("figures","sfig_simu_TypeI_omni.png"), 
       fig_simu_a_supp, width = 14, height = 5.6)

fig_simu_single_a_supp

# ggsave(file.path("figures","sfig_simu_single_omni.png"), fig_simu_single_a_supp, width = 10, height = 7)

fig_simu_b_supp

# ggsave(file.path("figures","sfig_simu_TypeI_mainpost.png"), fig_simu_b_supp, width = 13, height = 8)

fig_simu_single_b_supp

# ggsave(file.path("figures","sfig_simu_single_mainpost.png"), fig_simu_single_b_supp, width = 14, height = 14)

fig_simu_c_supp

Session information

In this project, we used R (Version 4.1.3; R Core Team 2022) and the R-packages afex (Version 1.1.1; Singmann et al. 2022), dplyr (Version 1.0.9; Wickham, François, et al. 2022), emmeans (Version 1.7.4.1; Lenth 2022), forcats (Version 0.5.1; Wickham 2021a), ggplot2 (Version 3.3.6; Wickham, Chang, et al. 2022), ggpubr (Version 0.4.0; Kassambara 2020), knitr (Version 1.39; Xie 2022), lme4 (Version 1.1.29; Bates et al. 2022), Matrix (Version 1.4.1; Bates, Maechler, and Jagan 2022), papaja (Version 0.1.1; Aust and Barth 2022), purrr (Version 0.3.4; Henry and Wickham 2020), RColorBrewer (Version 1.1.3; Neuwirth 2022), readr (Version 2.1.2; Wickham, Hester, and Bryan 2022), rmarkdown (Version 2.14; Allaire et al. 2022), stringr (Version 1.4.0; Wickham 2019), tibble (Version 3.1.7; Müller and Wickham 2022), tidyr (Version 1.2.0; Wickham and Girlich 2022), tidyverse (Version 1.3.1; Wickham 2021b), tinylabels (Version 0.2.3; Barth 2022), and vcd (Version 1.4.10; Meyer, Zeileis, and Hornik 2022).

sessionInfo()

## R version 4.1.3 (2022-03-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] emmeans_1.7.4-1    afex_1.1-1         lme4_1.1-29        Matrix_1.4-1      
##  [5] knitr_1.39         RColorBrewer_1.1-3 ggpubr_0.4.0       vcd_1.4-10        
##  [9] forcats_0.5.1      stringr_1.4.0      dplyr_1.0.9        purrr_0.3.4       
## [13] readr_2.1.2        tidyr_1.2.0        tibble_3.1.7       ggplot2_3.3.6     
## [17] tidyverse_1.3.1   
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-157        fs_1.5.2            lubridate_1.8.0    
##  [4] insight_0.17.1      httr_1.4.3          rprojroot_2.0.3    
##  [7] numDeriv_2016.8-1.1 tools_4.1.3         backports_1.4.1    
## [10] bslib_0.3.1         utf8_1.2.2          R6_2.5.1           
## [13] DBI_1.1.2           colorspace_2.0-3    withr_2.5.0        
## [16] tidyselect_1.1.2    rematch_1.0.1       compiler_4.1.3     
## [19] cli_3.3.0           rvest_1.0.2         xml2_1.3.3         
## [22] labeling_0.4.2      bookdown_0.26       bayestestR_0.12.1  
## [25] sass_0.4.1          scales_1.2.0        lmtest_0.9-40      
## [28] mvtnorm_1.1-3       digest_0.6.29       minqa_1.2.4        
## [31] rmarkdown_2.14      tinylabels_0.2.3    pkgconfig_2.0.3    
## [34] htmltools_0.5.2     highr_0.9           dbplyr_2.1.1       
## [37] fastmap_1.1.0       rlang_1.0.2         readxl_1.4.0       
## [40] rstudioapi_0.13     farver_2.1.0        ggh4x_0.2.1        
## [43] jquerylib_0.1.4     generics_0.1.2      zoo_1.8-10         
## [46] jsonlite_1.8.0      car_3.0-13          magrittr_2.0.3     
## [49] parameters_0.18.0   Rcpp_1.0.8.3        munsell_0.5.0      
## [52] fansi_1.0.3         abind_1.4-5         lifecycle_1.0.1    
## [55] stringi_1.7.6       yaml_2.3.5          carData_3.0-5      
## [58] MASS_7.3-57         plyr_1.8.7          parallel_4.1.3     
## [61] crayon_1.5.1        lattice_0.20-45     splines_4.1.3      
## [64] haven_2.5.0         cowplot_1.1.1       hms_1.1.1          
## [67] pillar_1.7.0        boot_1.3-28         estimability_1.3   
## [70] ggsignif_0.6.3      effectsize_0.7.0    reshape2_1.4.4     
## [73] reprex_2.0.1        glue_1.6.2          evaluate_0.15      
## [76] modelr_0.1.8        nloptr_2.0.0        vctrs_0.4.1        
## [79] tzdb_0.3.0          cellranger_1.1.0    gtable_0.3.0       
## [82] papaja_0.1.1        assertthat_0.2.1    datawizard_0.4.1   
## [85] xfun_0.31           xtable_1.8-4        broom_0.8.0        
## [88] coda_0.19-4         rstatix_0.7.0       lmerTest_3.1-3     
## [91] ellipsis_0.3.2      here_1.0.1

Reference

Allaire, JJ, Yihui Xie, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, Hadley Wickham, Joe Cheng, Winston Chang, and Richard Iannone. 2022. Rmarkdown: Dynamic Documents for r. https://CRAN.R-project.org/package=rmarkdown.

Aust, Frederik, and Marius Barth. 2022. Papaja: Prepare American Psychological Association Journal Articles with r Markdown. https://github.com/crsh/papaja.

Barth, Marius. 2022. Tinylabels: Lightweight Variable Labels. https://github.com/mariusbarth/tinylabels.

Bates, Douglas, Martin Maechler, Ben Bolker, and Steven Walker. 2022. Lme4: Linear Mixed-Effects Models Using Eigen and S4. https://github.com/lme4/lme4/.

Bates, Douglas, Martin Maechler, and Mikael Jagan. 2022. Matrix: Sparse and Dense Matrix Classes and Methods. https://CRAN.R-project.org/package=Matrix.

Henry, Lionel, and Hadley Wickham. 2020. Purrr: Functional Programming Tools. https://CRAN.R-project.org/package=purrr.

Kassambara, Alboukadel. 2020. Ggpubr: Ggplot2 Based Publication Ready Plots. https://rpkgs.datanovia.com/ggpubr/.

Lenth, Russell V. 2022. Emmeans: Estimated Marginal Means, Aka Least-Squares Means. https://github.com/rvlenth/emmeans.

Meyer, David, Achim Zeileis, and Kurt Hornik. 2022. Vcd: Visualizing Categorical Data. https://CRAN.R-project.org/package=vcd.

Müller, Kirill, and Hadley Wickham. 2022. Tibble: Simple Data Frames. https://CRAN.R-project.org/package=tibble.

Neuwirth, Erich. 2022. RColorBrewer: ColorBrewer Palettes. https://CRAN.R-project.org/package=RColorBrewer.

R Core Team. 2022. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Singmann, Henrik, Ben Bolker, Jake Westfall, Frederik Aust, and Mattan S. Ben-Shachar. 2022. Afex: Analysis of Factorial Experiments. https://CRAN.R-project.org/package=afex.

Wickham, Hadley. 2019. Stringr: Simple, Consistent Wrappers for Common String Operations. https://CRAN.R-project.org/package=stringr.

———. 2021a. Forcats: Tools for Working with Categorical Variables (Factors). https://CRAN.R-project.org/package=forcats.

———. 2021b. Tidyverse: Easily Install and Load the Tidyverse. https://CRAN.R-project.org/package=tidyverse.

Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, Hiroaki Yutani, and Dewey Dunnington. 2022. Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://CRAN.R-project.org/package=ggplot2.

Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2022. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.

Wickham, Hadley, and Maximilian Girlich. 2022. Tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.

Wickham, Hadley, Jim Hester, and Jennifer Bryan. 2022. Readr: Read Rectangular Text Data. https://CRAN.R-project.org/package=readr.

Xie, Yihui. 2022. Knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.